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00 ' We study the non-equilibrium statistical mechanics of a system of freely moving particles, in 

' which binary encounters lead either to an elastic collision or to the disappearance of the pair. Such 

, a system of ballistic annihilation therefore constantly looses particles. The dynamics of perturbations 

(N . 

around the free decay regime is investigated from the spectral properties of the linearized Boltzmann 
(~| ' operator, that characterize linear excitations on all time scales. The linearized Boltzmann equation 

, is solved in the hydrodynamic limit by a projection technique, which yields the evolution equations 

for the relevant coarse-grained fields and expressions for the transport coefficients. We finally present 
' the results of Molecular Dynamics simulations that validate the theoretical predictions. 
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I. INTRODUCTION 



Understanding the differences and similarities between a flow of macroscopic grains and that of an ordinary liquid 
I is an active field of research [1, 2]. From a fundamental perspective, it is tempting to draw a correspondence between 
the grains of the former and the atoms of the latter in order to make use of the powerful tools of statistical mechanics 
to derive a large scale description for the various fields of interest, such as the local density of grains. A key difference 
. between a granular system and an ordinary liquid is that collisions between macroscopic grains dissipate energy, due to 
' the redistribution of translational kinetic energy into internal modes. This simple fact has far reaching consequences 
I , [2, 3], but also poses an a priori serious problem concerning the validity of the procedure leading to the hydrodynamic 
' description. Indeed, the standard approach retains in the coarse-grained description only those fields associated with 
Ch ' quantities that are conserved in collisions (such as density and momentum). There is however good evidence -both 
numerical and theoretical- that in the granular case, a relevant description should include the kinetic temperature 
field, defined as the kinetic energy density ([2, 4] and references therein), which is therefore not associated to a 
conserved quantity. 

Our purpose here is to test a hydrodynamic description with suitable coarse-grained fields, for a model system where 
not only the kinetic energy is not conserved during binary encounters, but also the number of particles and the linear 
momentum. The ballistic annihilation model [5-10] provides a valuable candidate: in this model, each particle moves 
^vqj ' freely (ballistically) until it meets another particle; such binary encounters lead to the annihilation of the colliding 
' pair of particles. In addition, we introduce a parameter < p < 1 that may be thought of as a measure of the distance 
* \ to equilibrium, so that an ensemble of spherical particles in dimension d undergoing ballistic motion either annihilate 
■ upon contact (with probability p) or scatter elastically (with probability 1 — p). For the corresponding probabilistic 
' ballistic annihilation model, the Chapman-Enskog [11] scheme was applied recently [12]. The hydrodynamic equations 
I were derived and explicit formulas for the transport coefficients obtained. Our goal here is two-fold. First, we would 
^ ' like to shed light on the context and limitations of the derivation, by obtaining the hydrodynamic description directly 
. ^ from the linearized Boltzmann equation. Second, we aim at putting to the test the theoretical framework thereby 
obtained by careful comparison with numerical simulations of the annihilation process. For granular gas dynamics, 
\ the same program is quite complete, although challenges remain [1, 2]. The objective here is to initiate a similar 
. . . 1 formulation for the ballistic annihilation model in view of a more stringent test of the hydrodynamic machinery. 

The paper is organized as follows. We start in section II with a reminder of results derived in Refs [8, 9]. The 
kinetic description adopted is that of the Boltzmann equation, since it has been shown that for p = 1 (all collision 
events leading to annihilation), the underlying molecular chaos closure provides an exact description at long times, 
provided space dimension d is strictly larger than 1 [9]. We may assume that the same holds for an arbitrary but 
non vanishing value of p, since the density is then still a decreasing function of time. The focus is here on an 
unforced system, which is characterized by an algebraic decay with time of the total density and kinetic energy 
density (homogeneous decay state) [8, 9]. More precisely, we are interested in small perturbations around this state, 
so that the Boltzmann equation will be subsequently linearized. After having identified the operator that generates 
the dynamics of fluctuations, attention will be paid in section III to its spectral properties. This will provide the 
basis for flnding in section IV the evolution equations for the hydrodynamic flelds (i.e. those chosen for the coarse- 
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grained description) and for obtaining explicit formulas for the transport coefficients. Finally, our predictions will be 
confronted in section V against extensive Molecular Dynamics simulations. Such a comparison is an essential step in 
testing the foundations of the hydrodynamic treatment. 



II. THE BOLTZMANN EQUATION APPROACH TO THE HOMOGENEOUS DECAY STATE 

A. Non-linear description 

The Boltzmann equation describes the time evolution of the one particle distribution function /(r, Vi,i). For a 
system of smooth hard disks or spheres of mass m and diameter a, which annihilate with probability p or collide 
elastically with probability 1 — p, it has the form 

I + Vi • /(r, v„t)= pJa[f\f] + (1 - p)Jc[f\f], (1) 

where the annihilation operator is defined by [9] 

Ja[f\g]=-cr''-' Jdv2jd&e{v,2-&){^^i2-&)f{r,v,,t)g{r,V2,t). (2) 

The elastic collision operator Jc reads [13, 14] 

Jcim = a''-' Jdv2jd&e{v,2 ■ ^)(vi2 • &){b,' - l)/(r,Vi,i)5(r,V2,t), (3) 

with Vi2 = Vi — V2, 6 the Heavisidc step function, & a unit vector joining the centers of the two particles at contact 
and an operator replacing all the velocities vi and V2 appearing in its argument by their precollisional values 
and V2, given by 

6" Vi =vl = vi - (vi2 • o-)<T, (4) 

K^^2 = V2 = V2 + (Vi2 • &)&. (5) 

We assume that the system can be characterized macroscopically by coarse grained (hydrodynamic-like) fields, that 
we define as in standard Kinetic Theory in terms of the local velocity distribution function /(r, v, t) 

n{r,t) = jdvf{r,v,t), (6) 

n(r,t)u(r,i) = jdvvf{r,v,t), (7) 

^n{r,t)T{r,t) = J dv^V' f{r,y,t), (8) 

where n(r,t), u(r,t), and T(r,t) are the local number density, velocity, and temperature, respectively. We have 
introduced here V = v — u, the velocity of the particle relative to the local velocity flow. We stress that the 
temperature defined has a kinetic meaning only, but lacks a thermodynamic interpretation. It seems natural to 
consider these fields, as they arc the usual hydrodynaniical ficilds of the equilibrium system (with p = 0). It is however 
not obvious at this point that restricting our coarse-grained description to the above three fields provides a relevant 
and consistent framework. A major goal of this paper is to provide strong hints that this is indeed the case. We will 
show in particular that closed equations can be obtained for these fields in the appropriate time and length scales, 
under reasonable assumptions. 

The Boltzmann equation (1) admits a homogeneous scaling solution fn in which all the time dependence is em- 
bedded in the hydrodynamic fields, with the further simplification that those fields are position independent. The 
existence of this regime could not be shown rigorously, but, numerically, such a scaling solution quickly emerges from 
an arbitrary initial condition [8, 9]. It has the form [9] 

/H(v,i) = ^^Xh{c), with c = (9) 

VH{t) 
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where 

VH{t) =\''^X'' (10) 

is the "thermal" (root-mean-square) velocity and Xh{c) is an isotropic function depending only on the modulus c = |c| 
of the rescaled velocity. By taking moments in the Boltzmann equation and using the scaling (9), it can be seen that 
the homogeneous density and temperature obey the equations [12] 

dnH{t) ^ -p^^(t)CnnHit), (11) 



at 
dt 



-pVH{t)CTTH{t), (12) 



where we have introduced the collision frequency of the corresponding hard sphere fluid in equilibrium (with same 
temperature and density) 



d-l 



^hW- ^^T75 (d + 2)r(rf/2) ^^^> 

and the dimensionless decay rates C„ and (t, that are functionals of the distribution function 

PCn = ~ jdCijdC2T{ci,C2)XH{ci)XH{C2), (14) 
Kt = -| jdc,jdC2(^^-l^T{ci,C2)xH{Cl)XH{C2). (15) 

In these expressions, 7 is a quantity that does not depend on time, which reads 

_ 2nH{t)vH{t)a'^-' _ (d + 2)^/2^(d/2) 

and the binary collision operator T(ci,C2), that should not be confused with the temperature, takes the form 

T(ci,C2) = Jd&&{ci2-&){ci2-&)[{l-p)b-^ -I]. (17) 

Finally, we can write an equation for the scaled distribution function xh (c) in terms of the coefficients and operators 
defined above 



P 



d ' 

{dCr - 2Cn) + CtCi • -T— 

OCi 



X-h(ci)=7 dc2T{ci,C2)xHici)xHic2)- (18) 



The operator in the last equation is defined again by equation (4), but substituting (vi,V2) by (ci,C2). 

Although an exact and explicit solution of equation (18) is not known, its behavior at large and small velocities has 
been determined [8, 9]. In this work we will use the approximate form of the distribution function in the so-called 
first Sonine approximation (an expansion around a Gaussian functional form, see Appendix A), which is valid for 
velocities in the thermal region, and all the functionals of Xh{c), that is the decay rates and the transport coefficients, 
will be evaluated in this approximation [8, 15]. 

B. Linearized Boltzmann Equation 

In the remainder, we consider a situation where the system is very close to the homogeneous decay state, so that 
we can write 

/(r,vi,t) = /H(vi,i)+^/(r,vi,t), |5/(r,vi,t)|«/H(vi,i). (19) 
Substitution of equation (19) into the Boltzmann equation (1), keeping only linear terms in Sf, yields 

Ir+vi-v) Sf{r,v,,t) 



at 

PUamfn] + Ja[fH\Sf]} + {1-P) UcmfH] + Jc[/i/|<5/]} , (20) 
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Given the scaling form of / (Eq. (9)), it is convenient to introduce as well the scaled deviation of the distribution 
function, i^x, as follows 

5/(r,vi,t) = ^^(5x(r,ci,T) . (21) 
Moreover, Eqs. (11) and (12) suggest to use the dimensionless time scale r defined by 

T=\j\t'vH{t'), (22) 

which counts the number of collision per particle in the time interval [0,t]. Combining Eq. (13) together with Eq. 
(11, 12) yields immediately vuit) = {1/vh{^) + p{Cn + Qt / '2-)t)-'^ and thus 

^ = r^Or ^^t \ + Z^H(0)p(Cn + CT/2)t] • (23) 

P(2U + Ct) 

In this time scale r (22), these equations (11, 12) are easily integrated, yielding 

nuij) = (0) exp(-2K„T), (r) = T^(0) exp(-2KTT), (24) 

and power law behaviors in time ni/(t) oc t-2C„/(2C„+CT) a,nd Tnkp) oc t-2CT/(2C,i+CT) ^i^^ large time f^X . It proves 
also convenient to introduce Fourier components [with the notation /ik = /dr exp"*"^ "" ft,(r)] so that the evolution 
equation for a general k component of 5x is, in the r timescale, 

d 

— (5Xk(ci, r) = [A(ci) - iIh{t)\!. ■ ci] (5xk(ci, r). (25) 

In this equation, the time dependent length scale Ih = 2vh{t)/vh{t) is proportional to the instantaneous mean free 
path {Ih{t) (X. n^{T), see Eq. (13)) and the homogeneous scaled Boltzmann linear operator reads 



A(ci)/l(ci) = 7 yrfC2T(ci,C2)(l+Pl2)XH(cl)/l(C2) 



+p{2Cn - dCrMci) - pCtCi ■ ^h{ci). (26) 

CCi 

In this expression, the permutation operator P12 interchanges the labels of particles 1 and 2 and subsequently allows 
for more compact notations. In the present representation, all the time dependence due to the reference state is 
absorbed in the mean free path, obtained from Ih{t) oc n~^{T) as 

^h(t) =;ff(0)exp(2KnT), (27) 

which, as expected, is an increasing function of time. 

C. Linearized Hydrodynamic Equations around the homogeneous decay state 

Let us define the relative deviations of the hydrodynamic fields from their homogeneous values by 

p(r,r)^^^^ = /dc<5x(r,c,T), (28) 

nH{T) J 

^{v,t)=^-^^ = [dcc5x{v,c,T), (29) 

^^"■'^^^^1? = y'dc(^^-l)(5x(r,c,r), (30) 



where 5y{r,T) = T/(r, r) — yait) denotes the deviation of a local macroscopic variable, ?y(r,T), from its homogeneous 
decay state value, yH(t)- Taking velocity moments in the Boltzmann equation (25), we obtain the linearized balance 
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equation for the k components of the hydro dynamic fields 

d_ 

d_ 



-K2Cn + CT) 



Wk 



+ (T)k(pk + ^k) + ilH{T)k ■ n[(5xk] - P SCu[6xk] 



0, 



^-2p(C„ + Ct) 



^k 



-2KTPk + i^Zi/(T)k • (wk + </»[^Xk]) - P SCtISx^] = 0. 



Here, we have introduced the traceless pressure tensor and the heat flux as 

n[^Xk] = JdcAic)Sxuic,T), 



where A and S are defined as 



0[<5xk] = 

Ay(c) 
E(c) 



dcS(c)(5xk(c,r), 



2 



d+2 



and the functionals 



P 5Cn[5x] =7 jdCijdC2T{ci,C2){l+'Pi2)XH{ci)SXk{C2,T), 
P Ku[Sx] =7 JdCiJdC2CiT{ci,C2){l+Vl2)XH{ci)SXk{C2,T), 

P SCt[Sx]=1 jdcijdc2 - r(ci,C2)(l +Pi2)XH(ci)5xk(c2,r). 



(31) 

(32) 

(33) 

(34) 
(35) 

(36) 
(37) 

(38) 
(39) 
(40) 



The previous analysis therefore amounts to obtaining a set of complicated equations expressing the evolution of the 
hydrodynamic fields as a function of the rescaled homogeneous distribution function xm and the perturbation Sx- In 
order to obtain a closed set of equations for the hydrodynamic fields (31), (32), (33), we need therefore to express the 
functionals 11, 0, d(^n, ^Cu and S(t, in terms of the hydrodynamic fields themselves. We will see in the next section 
that, as long as we can treat Z/f (T)k as a small parameter and if the linear Boltzmann operator has some specific 
properties, it is possible to carry out this program and to close the linear hydrodynamic equations. However, since 
the mean free path ^//(t) increases with time (Eq. (27)), the requirement of a small Z//(T)k is necessarily limited to a 
time window depending on both k and the probability of annihilation p. An upper bound for this window is provided 
by the time when the mean free path becomes of the order of the system size. 



III. SOLUTION OF THE LINEARIZED BOLTZMANN EQUATION 

In this Section we explore the solutions to the linearized Boltzmann equation (25) and establish some properties 
of the homogeneous linear Boltzmann operator that will be essential for the coarse-grained description. Prom the 
expression of the linearized Boltzmann equation, we can identify the operator A — ik • cIh{t) as the "generator of the 
dynamic" of Jxk- As we are interested in the solutions of this equation in the hydrodynamic regime (large enough 
scales), it is convenient to study first the eigenvalue problem of the homogeneous linear Boltzmann operator. The 
inhomogeneous term will be treated perturbatively later on. 
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A. Hydrodynamic Eigenfunctions of A 



Let us consider the eigenvalue problem of the homogeneous linear Boltzmann operator A 



(41) 



Finding all the solutions of this equation is an insurmountable task. Nevertheless, it is possible to obtain some 
particular solutions, which will turn out to be the relevant ones in the hydrodynamic regime. The problem will be 
posed in a Hilbert space of functions of c with scalar product given by 

{g\h)= jdcxH\c)g*{c)h{c), (42) 

where g* denotes the complex conjugate of g. 

Of particular interest here are the eigenfunctions and eigenvalues associated with linear hydrodynamics. Following 
[16, 17], we use the fact that the homogeneous decay state is parameterized by the hydrodynamic fields uh, Th and 
\iH- Writing the Boltzmann equation satisfied by xh and differentiating it with respect to these fields allows then to 
obtain three exact relations from which one can extract eigenfunctions of the linearized Boltzmann collision operator. 
In Appendix B, we show that the functions 



d 

6(c) = Xh(c) + ^ • [cxh(c)] , 
d 

6(c) = zxh{c)- — -[cxHic)], 
d 

6(c) = -^Xh(c), 



with z = 2(n/CT, are solutions of Eq. (41), with eigenvalues 



Ai =0, 



Ao 



-p(Ct + 2C„), A3=Kt, 



(43) 
(44) 
(45) 

(46) 



respectively, A3 being d-fold degenerate. Although we cannot prove in general that these eigenvalues are indeed the 
hydrodynamic ones (i.e the upper part of the spectrum), we will assume that this is the case ; the self-consistency of 
the approach and comparison with numerical simulations will validate this assumption. Interestingly, in the particular 
case of Maxwell molecules where the full spectrum of A may be computed exactly (sec Appendix C), it appears that 
the above "hydrodynamic" modes dominate at long times, provided that p < 1/4. For larger values of p, the "kinetic" 
mode with largest eigenvalue decays slower than one of the three "hydrodynamic" modes. 

As a consequence of the non-hermitian character of the operator A, the functions {^/3}/3=i,...,3 are not orthogonal 
with respect to the scalar product defined in (42). They are nevertheless independent and, in order to define the 
projection onto the subspace spanned by these functions, it is necessary to introduce a set of functions {^/3}/3=i,...,3 
verifying the biorthonormality condition 



(?/3|?/3') = ^/3,/3'- 

Although the set {63}/3=i,...,3 is not unique, a convenient choice is given by 

2 + z z 



6(c) = 

6(c) = 

|3(C) = 



2{l + z) 
1 

2{l + z) 
cxh(c). 



1 + z d 



+ 



1 



1 + z d 



Xh{c), 
Xh{c), 



(47) 

(48) 

(49) 
(50) 



Indeed, the functions {63}/3=i,...,3 have to be linear combinations of xh{c), cxh{c) and c^xh{c} to ensure that 
projection of Sxk onto the {6} yields the coarse-grained fields pk, 6k and Wk, or combinations thereof. The fimctions 
{^/3}/3=i,...,3 span a dual subspace of that spanned by the eigenfunctions and for any linear combination of the 
hydrodynamic modes 



5(c) = ^a0^0{c), 

0=1 



(51) 
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the coefficients are given by 

af3 = {M = JdcxH\c)^0ic)g{c). (52) 

In particular, the projection of the distribution function Sxk on the subspace spanned by the functions is given by 
the coefficients 

mSx^)} = { Y^Pk - ^^f^^ Y^Pk + w,| . (53) 

Notably, these coefficients are simply linear combinations of the hydrodynamic fields linearized around the homoge- 
neous decay state. 

B. Projection of the Linearized Boltzmann Equation on the Hydrodynamic Subspace 

In this section, we study the linearized Boltzmann equation on the hydrodynamic subspace. Let us define the 

projectors 

3 

Ph{c) = J^i^pm^ic), (54) 

/3=1 



and 



so that any function can be decomposed as 



Pl = 1 - P. (55) 



h{c) = Ph{c) + P±h{c). (56) 



In the definition (54) we are considering the functions (43)-(45) and (48)-(50) defined above. 

Let us now consider the function Sxk- If we apply the projectors P and P± to equation (25), we obtain the following 
relations 



^ - P(A - tlH-k ■ c)P 



^-P^{K-ilH\i-c)Pi_ 



PSXi. = -PUh^ ■ cPx<5xk + P^P±Sxk, (57) 
P±Sxk = -Pi^Uh^ ■ cPSxu, (58) 



where we have used that 

P_iAP = , (59) 

which is obtained straightforwardly since are right-eigenfunctions of A. We note however that the {C/3}/3=i,...,3 are 
not left-eigenfunctions of A, so that PAPj_ 7^ 0. This also means that P and A do not commute. 

Equations (57) and (58) for the functions Pi5xk and P_Li5xk are coupled. Nevertheless, we shall see that, under 
certain conditions, we can decouple the equation for P(5xk in the long time limit. If we solve formally equation (58), 
we obtain 



P±'5xk(c,r) =Go(r)Pi(5xk(c,0)- / dr'Gr'ir - T')P±ilH{T')k ■ cP5x^{c,t'), 

Jo 



(60) 



where we have introduced the operator Gr'{T — r') defined from 

■^Gr'{T-T')=P^[A{c)-llH{T)k-c]P^Gr'{T-T'), G,,(0) = 1. (61) 

or 

If the hydrodynamic eigenvalues of the operator A are separated enough from the rest of the spectrum, the first term 
on the right hand side of (60) decays with the "non hydrodynamic" modes, faster than the second one. We can then 
write 

P±<5Xk(c,T) « - / dr'Gr'ir - T')P^ilH{T - T')k • cP5xk(c,T - r'). (62) 
Jo 
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and wc sec, by substituting equation (62) in (57), that we obtain an involved but closed equation for PSxk- It is 
worth pointing out that we have not proved scale separation, but assumed it in order to derive (62). For an explicit 
discussion of the scale separation assumption in a similar but somewhat simplified context, we refer to Appendix C, 
already alluded to above. 

The set of hydrodynamic equations (31)-(33) have been obtained through the projection of the Boltzmann equation 
onto the hydrodynamic subspacc. It now appears that, in the hydrodynamic time scale, the use of Eq. (62) will allow 
us to close these equations by substituting the distribution function by its decomposition in terms of the projectors, 
5xk = PSxk + -PL<^Xk- This is the aim of the next section. 



IV. LINEAR HYDRODYNAMIC EQUATIONS IN NAVIER-STOKES ORDER 

In this section we shall use the decomposition of Sxk into its hydrodynamic part, P6x\c, and non-hydrodynamic 
part, P±5xk, to close the linear hydrodynamic equations (31)-(33). We shall do so in Navier-Stokes order, that is, in 
the long time limit and in second order in the gradients (order fc^). 

Let us first introduce Pdxk in the linear pressure tensor and in the heat flux vector. Here the calculation is 
straightforward and we obtain 

U[PSxk] = 0, 0[P<5xk] = 0, (63) 

because the functions X-ff(c)A(c) and x/f(c)S(c) arc orthogonal to the subspacc spanned by the hydrodynamic 
eigenfunctions {^/3(c)}/3=i,...3. Turning our attention to the other functionals d(n, ^Cu and 6(t, the calculations 
become somewhat lengthy, and we show the details in the Appendix D. We obtain 

SUPSXk] = -^CnPk-CnOk, (64) 

SCu[PSxk] = -2C„Wk, (65) 
KT[PSxk] = -4CTPk- (3CT + 2C„)ek. (66) 

The negative signs occurring on the right hand side of these relations account for the fact that a fluctuation with a 
local enhanced density will induce an increased collision rate, hence a faster density decay. The same remark holds 
for temperature or local velocity flow fluctuations. 

We now have to calculate the contribution of Pi(5xk to the same functionals, to second order in k. This requires the 
knowledge of P^_5xk to first order in k since the heat flux and pressure tensor enter the balance equations (31)-(33) 
through their gradients and arc already weighted by a factor k. However, it should be noted that for consistency, the 
decay rates should be computed to second order in the gradients (see Eqs. (31)-(33)). We shall nevertheless restrict 
to first order, henceforth neglecting the various terms of order two that symmetry allows (such as V^n and V^T for 
5C,n and 5C,tj or as V^u for (^Cm)- We will further comment this approximation below. To leading order, we have that 
Gt-t'{t') « g^P±^P±T ^ so that we obtain from equation (62) 

-PL'^Xk(c, r) K, 

dT'e^^^^^'''P±ilH{T - T')k • cP^Xk(c, T - r') 

« -Ih{t) [ dr'e^^^-P^^'e-^P^-^'Pj.ik • cP(5xk(c, t - r'), (67) 
Jo 

where wc have used that Ih(t) oc e'^P'^"^ (Eq. (27)). We now have to relate PSxki^ ^ t') to PSxm{'^)j a-nd to be 
consistent with the approximation made above, we also have to restrict to leading order in k. In doing so, Markovian 
equations for the fields will be derived. From Eq. (57), we get 

3 

P5Xk(c,T-r')«e-^^^-'P5Xk(c,r) = ^e-^''-'(^^|JXk(T))C/3(c). (68) 

/3=1 

Substituting (68) in (67), we obtain an equation for P±Sxk to first order in k 

P±SxI'\c,t) = -Ih{t) Y.m6xk{T)) / dT'e^-(^-2fC"-^'')^--'p^zk. c$0(c), (69) 
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where /'I'^Xk = P^^Xu^ + 0{k'^). The pressure tensor and the heat flux up to first order in the gradients of the 
fields are now obtained by substituting equation (69) into equations (34) and (35). Taking into account the symmetry 
properties of the system, the resulting expressions can be written in the form 



0[^±^xL'^] = -i/H(T)k[K(r)^k+M(T)pk] 



(70) 
(71) 



Equation (70) is the expected Navier-Stokes expression for the pressure tensor, involving the shear viscosity coefficient 
fj, but equation (71) contains, besides the usual Fourier law characterized by the heat conductivity k, an additional 
contribution proportional to the density gradient and with an associated transport coefficient jl. This latter term is 
analogous to the one appearing in granular gases [3, 18]. 

The expression of the (time dependent) transport coefficients are 



flir) ^ JdcA^y{c)F3^^y{c,T) = ^'^^A^ (c)i^3,,j (c, r), 



(72) 



AW = ^^^^ ydcS(c) [Fi(c, r) + F2(c, r)] , (73) 
^^^^ = 2d{l + z) Z'^^^^^) [-^Fi(c,t) +F2(c,t)] , (74) 



where we have introduced the functions 



F3,ij{c,T) = /Ve^-(^-2fC''-fC-)-'p^c.6,,(c), 
Jo 



(75) 



Fi(c,r) 



/Ve^^(^-2f«")^'p^ca(c), (76) 

Jo 



F2(c,t) = / dr'e^^(^+^";^)"'PiC^2(c), 
Jo 



(77) 



and in the second equality of equation (72), we have summed over all the i, j, taking into account the symmetry of 
the linearized Boltzmann operator. 

Similarly, we calculate the deviations of the decay rates to first order in the gradients of the fields by substituting 
equation (69) into equations (38), (39) and (40). Taking into account the symmetry properties, we arrive at 

SUP±Sx'^^] = 0, (78) 
SUP^^X^^^] = ilH{TMCuA^)l^ + Cu,e{r)Ok], (79) 
SCT[P±dx'^^] = 0. (80) 

The expression for the coefficients are 

C«,p(t) = -^^Y^ Z"^^' JdC2XH{Cl)Ci2ici + C2) ■ [Fi (cj, t) + F2 (C2, r)] , (81) 



Cu,0{t) = 2^ JdClJdC2XH{ci)ci2{ci +C2)- [-2;Fi(c2,t) +F2(C2,t)] , 



(82) 



with P = 7r('^-i)/Vr[(d + l)/2] the d-dimensional solid angle. 

At this point, it is important to note that the transport coefficients defined in equations (72)-(74) and (81)-(82) 
are time-dependent, and this dependence is governed by the F^ functions. The (exponential) integrands appearing 
in the definitions of the Fj decay with the non-hydrodynamic (kinetic) modes, as a consequence of the action of 
the projector P±. Prom our assumption of hydrodynamic versus kinetic scale separation, all kinetic eigenvalues are 
smaller than the smallest hydrodynamic eigenvalue of A, which is A2 = —pCr — '2pCn- This ensures the convergence of 
the integrals (75)-(77) for r 00. In order for the transport coefficients to reach their r ^ 00 limit faster than any of 
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the hydrodynamic time scales, we need moreover the more stringent condition that the fastest kinetic mode is at least 
separated by a p(t gap from A2 = — pCt ~ 2pC„ : under this condition, the time dependence of the exponential term 
in the integral giving F2 is fast enough so that the transport coefficients, that depend on the functions through 
(72)- (74), can be considered as constants on the hydrodynamic time scale. With this proviso in mind, it is possible 
to set r ^ 00 in the integrals (75)-(77) and the time-independent transport coefficients obtained in this section are 
then equivalent to those calculated in reference [12] by the Chapman-Enskog method. We recall in Appendix A their 
expressions in the first order Sonine approximation. 

Finally, if we substitute the expressions derived above for the fluxes, equations (70)-(71), and the decay rates, 
equations (78)- (80), and we take into account that in the hydrodynamic time scale wc can substitute all the coefficients 
by their r ^ 00 limit, we obtain the following closed equations for the linear deviation of the hydrodynamic fields 



d \ 

— + 2pCn j PU + ilH{T)kWk\\ +pCnOu = 0, 



Wk_L = 0, 



d , 2(d- l),o , ,-,2' 



Wk|| 



+-iHiT)k [(1 - 2Ku,p)Pk + (1 - 2pCuM = 0, 



^+KT + ^/H(T);5fc2 



2Kt + -^kk^y^ik^ 



PVl^- i-^H{T)kw-^\\ = 0, 



where | and Wkj. are the longitudinal and transversal parts of the velocity vector defined by 

«'k||=Wk-k, Wk_L = Wk - Wkiik, 



(83) 

(84) 

(85) 
(86) 
(87) 



and k is the unit vector along the direction given by k. 

Equation (84) for the shear mode is decoupled from the other equations and can be readily integrated. If we 
introduce a non-dimensional wave number k = iij(0)fc, scaled by the mean free path at the time origin, we obtain the 
explicit solution 



Wk_L(r) = exp 



fjk^ 



Wk±(0). 



(88) 



Interestingly, depending on fc, the perturbation may initially increase if pCx ~ fjk'^ > 0. For long times however, the 
exponential e'*^^"'^ always dominates the linear term pCtt and the perturbation decays. 

The other three fields, namely, the density pk, temperature 9^, and the longitudinal velocity WkH, obey the system 
of coupled linear equations 





where the time-dependent matrix is 



M(r) 



-2pCn 

-^lH{T)kil-2pCu,p) PCt 
" " I,-, 1,2 



M(t) 



-ilH{T)k 



-2pCt - iiii{T)fik 



-i-^lH{T)k 



(89) 



■^lH{T)k{l-2pCu,e) 



2,2 



r.;^2 



lfj{T)kk 



(90) 



We note here that this matrix differs from Eq. (59) of Ref. [12], where the analysis amounts to overlooking the time 
dependence of the mean free path, so that all entries of the hydrodynamic matrix exhibit the same time dependence. 

The different time dependences present in Eq. (90) render the stability analysis more difficult. For long times however, 
the eigenvalues of the matrix M(t) are always negative and the perturbations a priori decay. A caveat is nevertheless 
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^ {n{0)/n{x)f 

Figure 1: A) Time evolution of a linear perturbation of the transversal velocity for a system with annihilation probability 
p = 0.1. The solid lines are the Molecular Dynamics results and the dashed lines are the theoretical predictions (no adjustable 
parameter). Note how the theoretical predictions correctly account for the increase of the perturbation at short times. The 
inset shows the evolution of density with rescaled time r, where a is the discs' radius. For r = 20, na^ ~ 9.10"*. B) 
Evolution of 5iikx(T)/(5ukx(0) as a function of {0) / (r) . The dashed line is the exponential decay predicted by Eq. (92). 
n^(0)/n^(r) = 400 corresponds to r ~ 15. 

in order. It is worth pointing out that equations (88) and (89) break down at long times. Once the function Z//(T)fc 
exceeds unity indeed, the expansion in the gradients we have performed is no longer vahd, and it would be necessary 
to include terms of higher order in k. In addition, if the perturbation initially increases sufficiently to leave the linear 
regime, our description breaks down and it becomes necessary to consider non linear terms. 

Although quite involved, the evolution equations (89) can be numerically integrated, using for instance the transport 
coefficients computed in the Sonine approximation in Appendix A. This is what we do in the next section, in order 
to compare the theoretical predictions with Molecular Dynamics simulations. 

V. MOLECULAR DYNAMICS SIMULATIONS 

To put our theoretical predictions to the test, we have performed Molecular Dynamics (MD) simulations of a system 
of N smooth hard disks {d = 2) which undergo ballistic flights punctuated by coUisional events : at each collision, the 
discs annihilate with probability p ; otherwise, they collide elastically. The particles are localized in a square box of 
size L with periodic boundary conditions. An event driven algorithm [19] has been used and the initial density has 
been chosen low enough to be always in the dilute limit. The parameters for all the MD simulations are A''(0) = 10^, 
reduced density n(0)cr^ = 0.05 where a is the discs' radius and < p < 1. The initial conditions we have considered 
correspond to small amplitude perturbations around the homogeneous decay state, to enforce the validity of the 
linearized hydrodynamic equations (83)-(86). 

A. Perturbation of the transversal velocity 

Since Eq. (84) for the shear mode is decoupled from the other equations, one of the simplest macroscopic pertur- 
bation one can think of consists in an initial harmonic perturbation of the transversal component of the velocity field, 
whose evolution is given by Eq. (88). We shall consider a small perturbation in real space of the form 

u^{r,0) ^ Asm{kmy), (91) 

with A — 1Q~^vh{0) and km — Svr/L, where L is the linear size of the system. The reason for choosing the smallest 
possible value of k compatible with the boundary conditions is twofold. First, the corresponding mode is the most 
unstable at short times (see Eq. (88)). For the parameters of the simulations, we indeed probe the region where 
pC,T > fjk^, so that the hydrodynamic equation (88) predicts an initial increase of the perturbation. Second, the low 
k regime is that where our large scale predictions are most likely to be relevant. 

Figure 1 displays the evolution of Wk^± as a function of the number of collisions per particle r for a system with 
p = 0.1, averaging data over 50 different trajectories. The solid line is the simulation result and the dashed line is 
the theoretical prediction, Eq. (88), where the shear viscosity and the decay rates have been computed using the 
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^ («(0)/«(x))' 

Figure 2: A) Same as in Fig. 1 but for a system with p = 0.5, which -loosely speaking- may therefore be considered as being 
more "distant" to equilibrium. The inset shows the evolution of density with rescaled time r. For r = 4, the rescaled density 
is ncr^ ~ 9.10"* (while n{Q)a^ = 0.05). B) Evolution of 5uk±(r)/(5ukx(0) as a function of n^(0)/n^(r). n^{Q)/n^{T) = 400 
corresponds to r ~ 3. 
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Figure 3: Symbols: ratio f]/C,n (normalized by its p ^ value) extracted from the exponential fit of 5uk±(T) to Eq. (92). The 
solid line is the theoretical prediction in the first Sonine approximation. 



standard tools of kinetic theory (here, the first Sonine approximation [12], see Appendix A). An excellent agreement 
is obtained without any adjustable parameter, including the predicted increase of the perturbation at short times, and 
as also observed for a larger annihilation probability p — Q.'b (see Fig. 2, in which the MD data are obtained by an 
average over 150 runs). 

Recalling that Jiif (T)/n/f (0) = exp(— 2pC„r) and that vh{t) — f^r(O) exp(— p^tt), it proves also convenient to 
consider the actual velocity field 5vl{t) = vh{t)'w{t) instead of its dimensionless counterpart w, since the prediction 
then takes the form : 



'5wk_L(r) exp ■ 



^^2 



T)k 



^iik±(0). 



(92) 



The plot of Sui^±{t) / Sui^±{0) as a function of (jt-h (0)/n//(T))^ allows then by simple exponential fitting to extract 
v/Cn (recall that k — ^^f(0)fc is known). Figure 3 compares the ratios fi/Cn extracted from such fits for various values 
of p with the theoretical prediction. We have plotted the value of the viscosity normalized by its elastic value, fje- For 
p = 0.1 the agreement is quite good but for p = 1 we obtain discrepancies of the order of 15%. Such deviations could 
be due to the limitation for high dissipation of the first Sonine approximation that has been used to compute the 
numerical values of the various quantities involved in the description (in particular, the deviations of the homogeneous 
decay state velocity distribution from its Gaussian form might be relevant). Additionally, the shear viscosity could 
suffer from finite size effects. For a related discussion in the realm of granular gases, where both effects alluded to 
are at work, see [20-22]. Finally, neglecting the contribution to SC,u might not be innocuous. From symmetry 
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Figure 4: Time evolution of the longitudinal velocity for a system with p = 0.1 as a function of r. The solid line shows the 
simulation results and the dashed line is the numerical solution of (89). For r = 15 the rescaled density is ncr ~ 2.5.10"^ The 
inset shows the increase of mean- free path with time and that for r > 10, Ink is no longer a small quantity. 
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Figure 5: Time evolution of ek„ and pk^ as a function of r for a system with p = 0.1. The solid lines show Molecular Dynamics 
results while the dashed lines are for the numerical solution of (89). For r = 15 the rescaled density is na^ ~ 2.5. 10~'^. 

considerations, such a term must be of the form fc^^_f/(T)^Wk, so that the equation (84) for the transversal vefocity 
has the same form as the one we considered, but with a "shifted" shear viscosity. It is worth pointing out here that in 
the corresponding equation (88), putative order k'^ corrections to S(n and S(^t play no role (see Eq. (84)) : the decay 
rates Cn and Ct appearing in (88) are fingerprints of the t dependence of Ih and of the rescaling procedure leading 
to Wk from the actual velocity flow. Those two decay rates are therefore properties of the homogeneous solution and 
do not suffer any finite k correction. 

B. Perturbation of the longitudinal velocity 

In order to investigate further the validity of the hydrodynamic equations, we consider a perturbation of the 
longitudinal velocity 

Ua;(r, 0) = Asm{kmx), (93) 

where A — 10^^11^(0) and k,n — 2n/L. Since the hydrodynamic matrix M depends on time and does not 
commute with M, we could not solve analytically the set of equations (89), and we turned to a numerical integration, 
using the transport coefficients computed in Appendix A. 

In Fig. 4, we have plotted the time evolution of the rescaled longitudinal velocity field, Wfe,„||, as a function of the 
internal time clock r (number of collisions per particle) for a system with p = 0.1. The results have been averaged over 
16 trajectories. It can be seen that the theoretical framework is able to account for the non trivial time dependence 
of the perturbation dynamics. Moreover, as the equations for 9k„^ and pk,„ are coupled with the equation for w^^n, a 
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Figure 6: Time evolution of the longitudinal velocity for a system with p — 0.5 as a function of r. The solid line shows the 
simulation results and the dashed line is for the numerical solution. We have also plotted with a dotted line the numerical 
solution considering the elastic values of k and /x (i.e. their limit when p — > 0"*"). The inset shows the increase of mean- free 
path with time. 
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Figure 7: Time evolution of and pkm as a function of r for a system with p = 0.5. The different lines have the same 
meaning as in Fig. 6. 

perturbation such as (93) induces a response of the above two other fields (at variance with the transversal velocity, 
whose dynamics is decoupled from the other three modes, at least at the linear level of description adopted here). In 
Fig. 5, we have plotted the energy density ek^ — dk^ + Pk^^ ^-nd pk^ as a function of r. The agreement with theory 
is very good both qualitatively and quantitatively, with an evolution that is well predicted till r ~ 15 ; for long times 
the simulation data become somewhat noisy (less statistics can be achieved due to the smaller number of particles 
left in the system. The number of particles at r = 15 is iV ~ 5000). 

In Fig.s 6 and 7, results for a system with p = 0.5 are shown. Data are averaged over 64 runs. In this case, 
the agreement is still good qualitatively, with similar shapes for the theoretical and numerical curves, but some 
discrepancies are observed. The most significant deviation from the theoretical prediction occurs for the density, 
Pk^ , where the value of the minimum and also its position are not predicted accurately. However, the hydrodynamic 
framework still captures correctly the trends of the complex dynamics of the perturbations. The above discrepancies 
could be ascribable to the failure of the first Sonine approximation for high dissipation (i.e. "high" p) or to finite size 
effects. We mention here that within the first Sonine approximation, the dimensionless coefficients p, and k exhibit a 
divergent behavior in the vicinity of p — 0.8 [12] which is presumably unphysical, and is an indication of the limitation 
of the method. For completeness and to assess the robustness of the predictions with respect to a modification of 
the numerical values of the key parameters, we have also reported in Fig.s 6 and 7 the predictions obtained when 
the transport coefficients take their elastic hard disc value (e.g. jl{p = 0.5) ~ 0.505 while p,{p = 0) = 0, as required 
by Fourier's law ). The important point is that the time evolutions are not significantly affected, and that the main 
features remain the same. Additionally, as mentioned after Eq. (66) and (92), a possible source of inaccuracy lies 
in the truncation of decay rates to their first order (fc^) in the gradients. While the corresponding terms have been 
shown to be small for inelastic hard spheres [23] (with the notable simplification there that the velocity decay rate 



15 



vanishes identically due to momentum conservation), their relevance in the present case has not been assessed, apart 
indirectly for the velocity decay rate (5(J„, by noting that second order corrections do not spoil the accuracy of the 
prediction (88), see Figures 1 and 2. 

VI. CONCLUSIONS 

The objective here has been to explore the validity of a hydrodynamic description based on density, momentum, 
and kinetic temperature fields for a gas composed of particles which annihilate with probability p or scatter elastically 
otherwise, by a direct analysis of the spectrum of the linear Boltzmann equation. The motivation mainly was to study 
the applicability of hydrodynamics to systems which a priori lack scale separation and in which there are no collisional 
invariants. The analysis performed here has been shown to lead to results equivalent to those obtained previously 
from the more formal Chapman- Enskog expansion [12], with the difference that the present approach considers linear 
excitations only. However, the current spectral method is arguably more straightforward and explicitly shows that 
the hydrodynamic description arises in the appropriate time scale, when the "kinetic modes" can be considered as 
negligible against the hydrodynamic ones. 

The eigenvalue problem of the Boltzmann operator linearized around the homogeneous decay state has been ad- 
dressed and we have identified the hydrodynamic eigenfunctions. These eigenfunctions are not simply linear com- 
binations of 1, V and as it happens in the elastic case, but they are replaced by derivatives of the homogeneous 
decay state velocity distribution function xh, that is not known analytically. As a consequence of the non-hermitian 
character of the linearized Boltzmann operator, the eigenfunctions are not orthogonal. It is nevertheless possible to 
construct a set of biorthonormal functions, {C/3}/3=i,...3, which are linear combinations of 1, v and v^, a crucial point in 
order to obtain the hydrodynamic equations. The analysis is complicated by the fact that none of the {^/g} fimctions 
are left eigenfunctions of the linearized Boltzmann operator, since no quantity is conserved during binary encounters. 
We have used these hydrodynamic eigenfunctions to derive to Navier-Stokes order the heat and momentum fluxes, 
together with the various decay rates. To this end, we have decomposed the distribution function, (5xkj into its 
hydrodynamic and non-hydrodynamic parts. This decomposition enables us to close the hydrodynamic equations in 
the long time limit and to order fc^, and provides Green-Kubo formulas for the transport coefficients. We then arrived 
at the linearized equations (around the homogeneous decay state) for the hydrodynamic fields, in the usual form of 
partial differential equations with coefficients that are independent of the space variable but depend on time, since 
the reference state considered is itself time dependent. If we analyze the stability of those equations, we may conclude 
that small perturbations should decay in the long time limit. Nevertheless, it must be stressed that the perturbation 
may increase at short times, thereby possibly leaving the linear domain where our analysis holds. The long time 
dynamics in such a case remains an open question. 

In section V, we have reported Molecular Dynamics simulations for the evolution of a perturbation of the transversal 
and longitudinal velocity fields, which show a rich dynamics. The agreement between theory and simulations is 
very good for moderate values of the annihilation probability p (say p < 0.5), which gives strong support to the 
theoretical analysis developed here. The theoretical curves still agree qualitatively at larger p, with however some 
quantitative discrepancies which might be a manifestation of the approximations underlying the computation of the 
transport coefficients (namely k and /i, evaluated to first order in a Sonine expansion [12]). Indeed, those coefficients 
are predicted to exhibit a divergent behavior for p ~ 0.8, but the simulations we have performed do not show a 
qualitatively different behavior for these values of p. A further complication -that is also an a priori limitation for 
the efficiency of a hydrodynamic approach- is that for values of p close to unity, the separation of time scales between 
the kinetic and hydrodynamic modes is not clear cut (the density decay rate is on the order of the collision frequency, 
comparable to the inverse typical time of the kinetic modes). To address this concern, one should study in detail the 
spectrum of non hydrodynamic modes to find the slowest, and compare it to the fastest decay rate in our problem (i.e. 
Cn). Such a program, left for future work, has been achieved in Appendix C within the Maxwell model framework, 
with the conclusion that scale separation does not hold for p > p* = 1/4. While the threshold p* is a priori model 
dependent, a similar phenomenon is to be expected in the original "hard-sphere" dynamics considered in this paper. 
The coarse-grained description for p> p* is an open question. 
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Appendix A: APPROXIMATE EXPRESSION FOR THE TRANSPORT COEFFICIENTS AND DECAY 

RATES 

For completeness, we recall here the explicit expressions used for the transport coefRcients and decay rates, as 
obtained within the first order Sonine scheme in Ref. [12]. The distribution function in this approximation reads 



d{d + 2) d+2 2 

8 2~^^ 



with the kurtosis 



a2 



8(3 - 2V2)p 



(4d + 6 - V2)p + 8V2(ci - 1)(1 - p) ' 
This expression allows us to calculate the decay rates 

. d+2f 1 



and also the transport coefBcients 



rf+2 / 8d+ll 



V = 



(Al) 
(A2) 

(A3) 
(A4) 

(A5) 



1 



1 , , d + 2,„ 
2KnM+ — ^(2a2 + l) 



(A6) 



2f„ - 3pCt - 2pCn 
where the values of the coefRcients Vn and are 

P 



KTK+^^(2a2 + l) 



,2 278 + 375d + 96(^2 + 2^3 
3 + 6d+2d"' -a2 — 



32(d + 2) 



+ (l-p) 1-02 



1 

'32 



(AT) 



(A8) 



32d 



16 + 27d + 8rf^ + 02 



,2 , 2880 + 1544rf - 2658^2 - 1539rf3 - 200d'* 



32d{d + 2) 



,d-l f 1 ' 



Finally, the expressions for (u,p and Cu.e are 



Cu, 



:-i) 



" (i(d+2) 



with 



{d + 2f 
32{d-l) 



1 + 02- 



rf(rf + 2) 

86 - lOld + 32^2 + 88rf3 + 28^^ 
32(d + 2) 



(A9) 

(AlO) 
(All) 

(A12) 
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Appendix B: EIGENFUNCTIONS OF A 

In this Appendix some of the details leading to the solution of the eigenvalue problem (41) are given. Consider first 
the function 

*i(c)=xh(c). (B1) 

and let the linearized operator A act onto xh 



A(ci)xh(ci) = 7 jdc2T{ci,C2){l + Pi2)xh{ci)xh{c2) 



+ p{2Cn - rfCT)XH(ci) -KtCi • ^^Xh{ci). (B2) 
If we take into account the equation for Xh{ci) 

d f 

{dCr - 2Cn)pXi?(ci) +KtCi • g^XHici) = 7 /rfC2T(ci,C2)XH(ci)XH(C2), (B3) 

we can rewrite equation (B2) as 

d 

A(ci)xi?(ci) = {dCr - Kn)pXH{ci) +KtCi • ^Xi?(ci). (B4) 



Consider now the function 



*2(c)=c-^Xh(c). (B5) 



In order to proceed, we perform the change of variables Ci = Xc'^ in equation (B3) 
(dCr - 2C„)pxh(Ac1)+Ktc'i-^Xh(Ac;) 



Deriving with respect to A we obtain 



d d f d 

{dCr - Kn)p-Q^XH {Xci) + KtCi • ^ \ dX^"^^'^^'> 



d 

{dCr - 2C„)p*2(ci) + KtCi • — *2(ci) 



(B6) 



= (rf+l)A'^7 jdC2T{ci,C2)XH{Xci)XH{XC2) 

+ 7A'^+^|dc2T(ci,C2)^^^^^XH(Ac2) 

+ jX'^+^ Jdc2T{cuC2)XH{Xc,)^^^^^. (B7) 
and taking A = 1 we arrive at the equation for ^'2(0) 



= (rf+l)7 jdC2T{ci,C2)XH{Xci)xH{XC2) 
+ 7 yciC2T(ci,C2)*2(Cl)X//(AC2) 

+ 7 ydC2T(ci,C2)XH(Aci)*2(C2), (B8) 
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or equivalently 



A(ci)f2(ci) = -(rf+l) 



d 

{dCr - 2C„)pXff(ci) +KtCi • g^XHici) 



(B9) 



Consequently, equations (B4) and (B9) can be written as 

A(ci)*i(ci) = (dCT-2C„)p*i(ci)+KT*2(ci), (BIO) 

A(C1)*2(C2) = -(d+l)[(rfCT-2C„)p*l(ci)+KT*2(ci)]. (Bll) 

With equations (BIO) and (Bll) we can easily sec that 

A(ci)[(d + 1)4-1(01) + *2(ci)] = 0, (B12) 

so that 

a(ci) = (d+l)*i(ci) + *2(ci), (B13) 

is an eigenfunction of the operator A(ci) with eigenvalue Ai = 0. It is also straightforward to see that 

6(ci) = {dCr - 2C„)?>*i(ci) +Kt*2(ci), (B14) 

is an eigenfunction of A. We obtain that 

A(ci)6(ci) = (dCT-2C„)pA(ci)*i(ci)+KTA(ci)*2(ci) 

= -(Ct + 2C„)p6(ci), (B15) 

where we have used equations (BIO) and (Bll). Therefore, we have that ^2(ci) is an eigenfunction of A(ci) with 
eigenvalue A2 = -(Ct + 2C„)p. 

Finally, let us consider the last function 

*3(c) = -^Xh(c). (B16) 

Deriving the equation obeyed by Xh{c — w), with respect to w and subsequently evaluating the result for w = 0, we 
obtain 

d 

(dCr-2C„)p*3(ci) + pCt*3(ci) +pCtCi • ^*3(ci) 
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ydC2T(ci,C2)(l + Pl2)*3(ci)Xff(C2), (B17) 



or equivalently 



A(ci)*3(ci) =Kt*3(ci). (B18) 
In other words, ^3(01) = *3(ci) is an eigenfunction of A(ci) with eigenvalue A3 = p(t- 

Appendix C: LINEARIZED BOLTZMANN OPERATOR FOR MAXWELL MOLECULES 

The objective in this Appendix is to study the spectrum of the linearized Boltzmann operator for Maxwell molecules 
with annihilation. It will be shown that for < p < p* , where p* depends on the specific Maxwell model under 
consideration, the norm of the hydrodynamic eigenvalues are smaller than the rest of the spectrum. 

The main characteristic of Maxwell models is that the difi'erential cross section multiplied by the relative velocity is 
independent of the relative velocity. We are going to assume that it is also independent of the angle between the two 
colliding particles. Then, the Boltzmann equation for a system of Maxwell molecules which annihilate in a collision 
with probability p and collide clastically otherwise (with probability 1 — p) reads 

(^+vi-v)/(r,vi,t) = -pl3n Jdv2f{r,Vi,t)f{r,V2,t) 

+ {l-p)f3 jdw2jda[b-^ -l\f{v,wut)f{v,w2,t), 

(CI) 
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where /? is a constant representing the microscopic scattering cohision frequency, fl = /T{d/2) is the d- 

dimensional sohd angle and the operator b'^^ is defined in the main text, Eq. (4). 

If we consider the homogeneous case, it is straightforward to see that the temperature is a constant in time, 
TH{t) = Th{0), and that the density evolves as 

1 + PCnt 

with Cn = l3flnH{0). Moreover, it was shown in [24] that there is an exact mapping between the homogeneous equation 
for Maxwell molecules with annihilation (arbitrary p) and the usual elastic Maxwell molecules {p = 0). If we introduce 
the time scale 



it can be seen that the distribution function for arbitrary p is 

f{^,t) = ^j''[^,i^-p)m, (C4) 

nH{Oj 

where nnit) is given by formula (C2), s{t) by (C3), and the hmction is the distribution function for elastic Max;well 
molecules. Note that this relation is also valid for p = 1 where is frozen in the initial condition. 

In the elastic case, every homogeneous distribution tends to relax to a Maxwellian after a transient time. Then, 
the same is going to happen for arbitrary p < 1 because of the mapping. In this sense, we can consider that the state 
analogous to the homogeneous decay state introduced in the main text for hard particles and that will constitute the 
appropriate reference will be characterized by the distribution function 

fH{v,t) = ^XM{v/VH), (C5) 

where Vh = (^^)^^^ and Xm(i') = l/7r''/^e~"^ is the Maxwellian distribution. Note that in the homogeneous decay 
state, both the density and temperature decay, whether in the present case only the density decays. Then, as vh 
plays no role, we will consider units with = 1 for simplicity. 

Let us study now the linear response to an inhomogeneous small perturbation around the reference state as it was 
done in section II. If we introduce the scaled distribution function 

Sx{r, V, r) = ^ [/(r, v, t) - /^,(v, t)], (C6) 
nH{t) 



^ + /iff(s)vi-V 



the equation for Sx in the s scale defined in (C3) is 

6x{r, V, s) = A(vi)(5x(r, v, s), (C7) 
where we have introduced the function /ih(s) = nH{0)/nH{s), and the linearized Boltzmann operator 

A(vi)5(vi) = jd^f2jda[h-^ -l]{l + Vi2)xM{vi)g{^2) 

- PCnXM{vi) jdV2g{V2). (C8) 

Let us stress that, although there is an exact mapping for the full non-linear homogeneous equation between p = Q 
and arbitrary p, no such mapping exists for the linear inhomogeneous Boltzmann equation, Eq. (C7). Then, as in 
the main text, the possibility of an hydrodynamic description depends on the properties of the linearized Boltzmann 
operator. Here we will see that it is possible to calculate all the eigenf unctions and eigenvalues of this operator and 
that there is a region of the parameter p in which we have an appropriate scale separation. 
Let us write the linearized Boltzmann operator as 

A(vi)ff(vi) = (1 -p)A^(vi)ff(vi) -KnXM(^;i) JdV2g{V2), (C9) 
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Figure 8: Hydrodynamic eigenvalues, Aoo, Aio, and Aoi, and the slowest kinetic eigenvalue, Afc, as a function of the dissipation 
parameter p. The eigenvalues are normalized by Cn = /3Onif(0). 



where we have introduced the hnearized Boltzrnann operator for elastic Maxwell particles 

A''(vi) = ^ Jdv^Jd&ib-' - l]{l+ri2)XM{vi)9{v2), (CIO) 
whose spectral properties arc well known [25, 26]. In particular, for rf = 3 its cigcnfunctions are 

'/'HmW = ArlXMiv)Sl^-^/2{'"^)v^Yi^{e, if), r = 0, 1, . . . (Cll) 

where yjm(0, (p) are the spherical harmonics, functions of the polar angles (6,(fi) of v with respect to an arbitrary 
direction, S^j^^^^i'^'^) ^'^^ the Sonine polynomials which satisfy 



' dxe--S^{x)S^'{x) = ^i!i±f±ll^„„,, (C12) 



and Ari are some constants that are introduced in order to normalize the eigenfunctions and that play no role in 
the following analysis. The eigenvalues of A^, A^, are also known. It can be seen that A^, A{^, and A^^ (which is 
3 times degenerate) vanish, corresponding to the five hydrodynamic eigenvalues and that the slowest kinetic mode 
corresponds to an eigenvalue A^" = — Cn/3- 

The important point here is that the functions <f)rim are also eigenfunctions of the second operator in (C9). Taking 
into account the orthogonality properties of the spherical harmonics and of the Sonine polynomials, Eq. (C12), one 
has 



rfv^r;m(v) = 5ro5lo. (CIS) 



Then, as 0ooo(v) = Xm{v), the functions (t>rim{^) are in fact the eigenfunctions of the total linearized Boltzmann 
operator for arbitrary p. With the aid of (C13) is straightforward to see that the eigenvalues are 

Ari = (1 - p)\fi - pCnSroSlO. (C14) 

In Fig. 8 we plot the hydrodynamic eigenvalues as well as the slowest kinetic eigenvalue, A^, as a function of the 
dissipation parameter p. It can be seen that for < p < 1 /4, there is scale separation in the sense that the three 
modes (density, linear momentum and kinetic energy) retained in the coarse-grained description decay slower that any 
of the other "kinetic" modes. On the other hand, for 1/4 < p < 1, the largest kinetic eigenvalue is slower than Aqo- 
Wc therefore conclude here that a conservative requirement for the validity of our approach in the case of Maxwell 
molecules would he p < 1/4. 
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Appendix D: EVALUATION OF 5^, dCu AND 5Ct 

In this Appendix we calculate the contribution of the hydrodynamic part of ^Xk to the functional S(n, ^Cu and 
5(t defined in (38)- (40). To this end, we write explicitly P6xk as 



•P^Xk(ci) = Pk6(ci 



2^k + Pk ) 6(ci) + Wk • 6(ci) 



1 ^ d 
= PkXH(ci) - 2^k^ • [ciXh(ci)] - Wk • q^Xh{ci). 

Let us first evaluate S(n[PSxk]- After some algebra it can be seen that 

SCnlxnici)] = -4C„, 



SCn 



d 

dci 



[ciXh(ci) 



d 

dci 



Xh{ci 



0, 



(Dl) 

(D2) 
(D3) 

(D4) 



where we have used equations (B5), (B8), the definition of the density decay rate, equation (14), and symmetry 
considerations. Then, if we consider equations (D2)-(D4) wc finally obtain 



(D5) 



Now let us calculate SCu[PSx^]- Using the definition of the density decay rate, equation (14), and symmetry 
considerations, it appears that 



SUXH (ci)] = 0, 



SCu 



dci 



(ciX//(ci)) 



d 



dc 



■Xi/(ci) 



= 0, 

= SijSCu 



d 
dcv 



■Xi/(ci) 



= 2C„. 



(D6) 
(D7) 

(D8) 



Then, we have 

SCu[P6x).] = -2CnWk. (D9) 
Finally, we turn to JCrf-P^Xk]- Using the definitions of the decay rates, equations (14) and (15), we obtain 

^Ct[xh(ci)] = -4Ct, (DIO) 

r r) 

SCt 



6(7 



;h— Xh(ci) 
aci 



6Ct + 4C„, 




(Dll) 
(D12) 



from which it follows that 



Kt[P5x^.\ = -4CTPk - (3Ct + Kn)Ow 



(D13) 
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